function dRdv = gradR(H,impossible,K,Nu,phi,rho,tau,P)
% Gradient of conditional value function with respect to 22 flow utility 
% variables. Returns a 22x1 gradient for each type and combination of 
% observed variables profile and each preference point

dRdv=cell(1,tau);
for idx2 = 1:tau
    
	dRdv{idx2} = zeros(H,K,Nu);
	
    for idx1 = 1:Nu
        
        P0 = P{idx2}(1:end-1,:,idx1);
		temp = P0;
        temp(impossible(1:end-1,:)) = 0;
		P0(2:end,:)=temp(2:end,:);
		
        phiP0 = phi(1:end-1,:).*P0;
		omphiP0 = (1 - phi(1:end-1,:)).*P0;
        P23 = P{idx2}(end,:,idx1);
		
		for idx3 = 1:H
			
            A = zeros(K);
            
			for idx4 = 1:K
            
                A(idx4,1) = rho*sum(phiP0(:,idx4))/(1 - rho*P23(idx4));
                A(idx4,min(idx4+1,K)) = rho*sum(omphiP0(:,idx4))/(1 - rho*P23(idx4));
                
            end
            
            B = (phiP0(idx3+1,:)./(1 - rho*P23))';
			dRdv{idx2}(idx3,:,idx1) = (eye(K) - A)\B;
            
		end		
    end
       
end

